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Abstract. Metabolic pathways describe chains of enzymatic reactions. 
Their modelling is a key point to understand living systems. An en- 
zymatic reaction is an interaction between one or several metabolites 
(substrates) and an enzyme (simple protein or enzymatic complex build 
of several subunits). In our Mitochondria in Silico Project, MitoScop, we 
study the metabolism of the mitochondria, an intra-cellular organelle. 
Many ordinary differential equation models are available in the literature. 
They well fit experimental results on flux values inside the metabolic 
pathways, but many parameters are difficult to transcribe with such 
models: localization of enzymes, rules about the reactions scheduler, etc 
Moreover, a model of a significant part of mitochondrial metabolism 
could become very complex and contain more than 50 equations. 
In this context, the multi-agents systems appear as an alternative to 
model the metabolic pathways. Firstly, we have looked after membrane 
design. The mitochondria is a particular case because the inner mito- 
chondrial space, ie matricial space, is delimited by two membranes: the 
inner and the outer one. In addition to matricial enzymes, other enzymes 
are located inside the membranes or in the inter-membrane space. Anal- 
ysis of mitochondrial metabolism must take into account this kind of 
architecture. 



1 Introduction 

Metabolic pathways are a set of reactions catalyzed by enzymes which takes 
place inside cells. They constitute the cellular metabolism. In biology, numer- 
ical simulations allow interactions between molecules to test hypotheses about 
normal or pathologic behaviours. Metabolism can be viewed as a network. We 
have studied mitochondria, a specific intra-cellular organelle, which has its own 
enzymes and thus its own network, generally analyzed independently of the cell 
metabolism. 

A classical way of modeling reaction chains is to use ordinary differential 
equations (ODE). With such models it is possible to describe the evolution of 



the system. But many biological aspects are difficult to describe with ODE: par- 
ticular situation where there are few molecules reacting (no average behaviour), 
influence of the mass or the conformation of the molecules, transient association 
of several molecules. However such examples correspond to biological cases like 
disturbance of the network caused by lack of some molecules during a small por- 
tion of time or genetic mutations which change enzyme conformation. All of this 
shows that new models for simulation like these using multi-agents systems 
may improve the understanding of biological processes. 

One way to design a biological network is to see it as a world filled with more 
or less complex entities. One entity could be a simple molecule like a metabolite 
or an enzymatic complex. Each reaction is described by interactions between 
entities belonging to this world. In our model, each entity is represented by 
a reactive agent. This agent is situated in a 3D space and the reactions are 
scheduled by time steps. 

In this paper, we show what is needed to model metabolism pathways with 
a multi-agents system (MAS) in the context of the mitochondrial metabolism. 
We present a first approach to a 3D molecular design and the application to the 
particular case of the two mitochondrial membranes: the inner and the outer 
one. Each of them is composed of a phospholipids bilayer containing enzymatic 
complexes. The attention is focused on a generic description of molecular inter- 
actions in order to use this behaviour both for simple case as lipids and for more 
complex one like respiratory chain enzymes. 



2 Mitochondria: their membranes and their energetic 
metabolism 

Mitochondria (see Fig.l) are intra-cellular organelles which are the power plant of 
the cell. They are delimited by a double membrane, composed of phospholipids, 
which are naturally arranged in a bilayer to minimize their interactions with 
water. The energy delivered by mitochondria comes from a complex mechanism, 
involving a series of oxido-reduction reactions called "Respiratory chain" using 
oxygen we breathe in and the nutriments we ingest. They are catalyzed by 
four highly structured protein complexes and electron transporters (CoenzymcQ 
and CytochromeC) embedded in the inner mitochondrial membrane. Electrons 
transfer is linked to an extrusion of protons outside the mitochondria, creating an 
electro-chemical gradient, which will be used for ATP synthesis and other works 
[6] [7]. This process involves the respiratory chain and ATP synthesis, both of 
them called " Oxidative phosphorylation" . Oxidative phosphorylation has been 
modelled for a long time in order to integrate all the kinetic and thermodynamic 
properties of chemiosmosis theory developed by Peter Mitchell [7] . 
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Fig. 1. Mitochondria from http:/ /universe-review. ca/ 

However some aspects, particularly the electron transports between com- 
plexes I or II and complex III of the respiratory chain by a molecule called 
CoQ presenting different redox states, is rather complex [8]. CoQ is certainly a 
central point in electron transport in respiratory chain. This transport appears 
rather difficult to model with ODE [I] due to the great variety of interactions, 
of the reactions involved at different reacting sites and of the amount of the var- 
ious redox CoQ species. In addition, these reactions take place inside the inner 
mitochondrial membrane which is a highly heterogeneous medium (see Fig. I). 

All these questions could be addressed with a multi-agents model. In fact these 
problems illustrate how an MAS approach is highly suited to take care of spatial 
disposition of the various enzymatic complexes, competition between different 
complexes to access to some finite resources, complex interactions between sev- 
eral types of molecules (membranes and enzymes) sometimes in low quantity, 
environment effects (locally or not). 

In the next section, we firstly describe the generic characteristics of the re- 
tained MAS. Finally, we describe the specific molecular modeling (section 4). 

3 A multi-agents model for metabolic pathways 

Our MAS uses a 3D continuous space where the agents are situated. Situated 
Agents use space as modality of their interactions. The distances between agents 
are used to describe neighbouring rules. All the agents have a local perception of 
their neighbourhood. A grid divides space into neighbouring units which simplify 
interactions complexity giving directly neighbour positions. 

The MAS is called open because of the variable number of the agent during 
simulation. The continuous degradation and synthesis of molecules is a common 
biological process. 



The time is discretized and a scheduler is defined for each type of agents. 
Actually, the inner-mitochondrial membrane is composed of enzymes and phos- 
pholipids. An enzymatic reaction should not be simulated on the same time-scale 
as two phospholipids interaction. One can also notice that it is possible to sim- 
ulate time continuity with a discretization based on events and not on time 
steps. 

In general, the level used for representing the biological objects, i.e. the model 
granularity, varies from atom to molecular groups (see [4] for a short review). 
One can notice that an abstraction by a single point is not sufficient to take into 
account a dynamic 3D structure and its spacial orientation (see Fig. 2). On the 
contrary, molecular modeling simulation (atom by atom) is not possible for a 
whole organelle due to the large number of entities to compute. 
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Fig. 2. Potentials deformations applied to a simple molecule. 



A phospholipids model published by Edwards [2] uses two interacting points 
per molecule. In order to design a more generic model, we have extended this pre- 
vious model. Our model is suitable for modeling metabolites and enzymes alike. 
Thus we have chosen to reduce each agent (molecule) to its gravity center and 
its interacting points. An interacting point is a portion of molecule that could be 
affected by external forces: hydrophobia, electrostatic charges, etc Fig. 3 shows 
the generic 3D structure (part A) and its application to a phospholipid design 
(part B). Forces and torques are applied at the gravity center. Forces induce 
linear acceleration whereas torques induce angular acceleration as in classical 
mechanics. 

Spatial localization of each interacting point is defined by a vector in the local 
frame associated with the agent (molecule) using its gravity center as origin. This 
design gives information about the 3D structure of the molecule and its space 
orientation but also its internal dynamics. 
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Fig. 3. A: abstracted 3D structure. B: application to a phospholipid. 



Modeling internal dynamics could be important parameters in enzymatic 
reactions. For example, an already filled active site will not be able to react with 
another metabolite. A lot of enzymes interact with more than one metabolite, 
in a certain order or randomly, before transforming them into new products. 
The enzyme conformation could evolve after the binding of a first substrate, 
revealing a second active site for a second type of substrate. Our model allows 
these flexibles movements. 

Even if their freedom is limited, the interacting points can be seen as agents 
and thus molecules as communities of linked agents. 

Each pair of interacting points defines one direction with two types of in- 
teraction - attraction/repulsion. Therefore, n interacting points could induce 
2n(n — 1) forces with intensity depending on the distance between agents. As 
shown in Fig. 4, the interactions between phospholipids can be given by 4 forces 
fii faifii /s that describe hydrophobia, another one fa for hydrophilic and elec- 
trostatic charges, and the two last ones fa, fa for the molecular incompressibility. 
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Fig. 4. Example of interaction forces for phospholipids. 



As a first step, we have implemented simple functions with three parameters 
(Fig. 5 A): etfc - maximal magnitude value of fk, - the radius value under which 
fk magnitude is maximal - and r' k - the radius value over which fk is null. A sec- 
ond step will be to study the impact of the choice of functions over phospholipid 
structure formation that would emerge from MAS simulations. Fig. 5 B shows 
another function possibility. 
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Fig. 5. Forces intensities are function of the distance. A: simple ramp function. 
B: example of other function available. 



With this approximation, the interactions between agents are reduced to a 
set of forces. For a molecule i in interaction with its neighbors, the compound 
force is: 



Fi(t) = a Hn VVV / h .k(t)+ b Un randi m (t) - c Un V l (t) 



v , 0lin Tiin 

where: 

> au n is the forces from the neighborhood (forces from the k interacting points 
of the molecule j acting upon the h interacting points of the molecule i), 

> Pun is the thermal effects depending on molecule type, 

> Jim is the friction (V* = molecule velocity) depending on molecule type. 

Linear and rotational accelerations are thus given by Newton's Second Law 
- the fundamental law of dynamics: the acceleration of an object of constant 
mass is proportional to the resultant force acting upon it ( F i(t) = rrfa). This 
compound force is used in the power series (order 2) of the position function 

%{t): 

%{t + st) = %(t) + ^T-st + ^|r^ 2 + o(st 2 ) 

AXi(t + St) = V l (t)St + ^~a l (t)St 2 + o(St 2 ) 
<^> AX l (t + St) = V t (t)St + ^-~Fi(t)5t 2 + o(St 2 ) 
Each time step of the simulation requires for the linear movement: 



AVi(t + 6t)K>±Fi{t)5t 

AXi(t + St) » Vi(t)5t + ^Fi{t)St 2 



where: 



> m is the molecule mass, 

> X{ is the vector position of the gravity center of the molecule (i.e.: the origin 
of the local frame). 

The forces applied by neighboring molecules induce torques. The compound 
torque is: 

Ti(t) = a rot V" V] ~Th.k(t) + b rot rand rot (t) - c rot Wi{t) 

^ V ^ ^ v J 

h jj^i k „ 



where: 



> a rot is the torques from the neighborhood (torques from the k interacting 
points of the molecule j acting upon the h interacting points of the molecule 

i). 

> f3 ro t is the thermal effects depending on molecule type, 

> 7 ro t is the rotational friction depending on molecule type. 

In the same way, each time step of the simulation requires for the angular 
movement: 

AWi(t + 6t) « jT t {t)5t 
A0i(t + St) w Wi(t)St + ±Ti(t)5t 2 

where: 

> I is the rotational inertia linked to the molecule size, 

> &i is the vector perpendicular to the rotation plan, oriented so that the 
rotational movement is positive and which has a norm corresponding to the 
angular variation. 

This model with interaction points and gravity center seems to be a right 
way to take into account the internal dynamics, the 3D structure and the space 
orientation of the different biological molecules. Moreover the possibility to adapt 
the agent granularity gives the possibility to design heterogeneous molecules. 

4 Application to the mitochondrial membranes and 
mitochondrial metabolism 

In a first attempt, we used this model to design mitochondrial membranes. Two 
embedded membranes appear as an important characteristic of mitochondria. 
The inner membrane is the more irregular one. As can be seen in Fig.l one 
can understand that giving a mathematical definition of this form will be rather 
difficult. On the other hand, using an emergent property of MAS (analogous to 
natural membrane formation) can presumably be a better solution. It can take 
into account the diversity of phospholipids constituting the bilayer and the fact 
that the inner mitochondrial membrane is also composed of a large amount of 
proteins complexes. 

MAS also appear well suited to simulate the respiratory chain reactions, par- 
ticularly the electron transport, which occurs inside the inner membrane and the 
proton transport which occurs across it. In principle, the mitochondrial applica- 
tion of our MAS model should allow the natural emergence of the complex elec- 
tron transport by CoenzymeQ (CoQ) described under the name of " Q-cycle" [8] , 
based on the reactions rules described for each CoQ site on complexIII. 

There is another field, in which MAS would certainly have a great interest ; it 
is the field of mitochondrial diseases due to deficiencies in respiratory complexes. 
Mitochondrial diseases are a group of complex pathologies, which can affect 



different unrelated tissues with an unexpected association of symptoms ([5], 
[10]). 

The disturbance in the distribution of the various CoQ species due to mu- 
tations in respiratory chain complexes (giving sometimes only mild respiratory 
chain deficiency) seems to be at the origin of (severe) mitochondrial diseases 
[9]. Our MAS model easily describes the distribution of the various redox CoQ 
species during the evolution of the reactions, for any state of the respiratory 
complexes. It thus will give a strong help to describe the plausible properties 
of the system and to design new experiments allowing us to better understand 
these complex pathologies. 

Finally reduced CoQ is an antioxidants [3]. Chemical derivatives of this 
molecule are studied in order to reinforce its antioxidant properties and to treat 
the oxidative stress in some mitochondrial diseases or during aging. Our MAS 
model will help specifying the relevant targets of these molecules, to keep them 
in their antioxidant form and thus help designing their properties. 
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